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The list of putative sources of gravitational waves possibly detected by the ongoing worldwide 
network of large scale interferometers has been continuously growing in the last years. For some of 
them, the detection is made difficult by the lack of a complete information about the expected signal. 
We concentrate on the case where the expected GW is a quasi-periodic frequency modulated signal 
i.e., a chirp. In this article, we address the question of detecting an a priori unknown GW chirp. 
We introduce a general chirp model and claim that it includes all physically realistic GW chirps. 
We produce a finite grid of tempfate waveforms which samples the resulting set of possible chirps. 
If we follow the classical approach (used for the detection of inspiralling binary chirps, for instance), 
we would build a bank of quadrature matched filters comparing the data to each of the templates of 
this grid. The detection would then be achieved by thresholding the output, the maximum giving 
the individual which best fits the data. In the present case, this exhaustive search is not tractabfe 
because of the very large number of templates in the grid. We show that the exhaustive search can be 
reformulated (using approximations) as a pattern search in the time-frequency plane. This motivates 
an approximate but feasible alternative solution which is clearly linked to the optimal one. The time- 
frequency representation and pattern search algorithm are fully determined by the reformulation. 
This contrasts with the other time-frequency based methods presented in the literature for the same 
problem, where these choices are justified by “ad-hoc” arguments. In particular, the time-frequency 
representation has to be unitary. Finally, we assess the performance, robustness and computational 
cost of the proposed method with several benchmarks using simulated data. 

PACS numbers: 04.80.Nn, 07.05.Kf, 95.55.Ym 


I. INTRODUCTION 

The worldwide network Q of large scale interferomet¬ 
ric gravitational wave (GW) detectors have started to 
take data. The network includes the detectors GEO600, 
LIGO and TAMA. It will be completed soon by the up¬ 
coming Virgo. The overall sensitivity of these detectors is 
continuously improving. Interesting upper-limits for the 
amplitude of GWs are being set and the first detection 
is hopefully not too far. 

A large variety of astrophysical sources are expected 
to emit GWs in the observational frequency bandwidth 
of these detectors. From the data analysis viewpoint, the 
detection methodology for these sources depends on the 
availability of a reliable and complete model of the GW. 

Generally speaking, the oscillations of the GWs are 
related to the orbital, rotational bulk motion of the con¬ 
stituents of the emitting system. Since the system loses 
energy by radiation, or because of some other physical 
process involved, its orbital period, and consequently the 
GW frequency can vary with time. In such case, the emit- 
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ted GW is a frequency modulated signal i.e., a chirp. A 
detailed knowledge of the dynamics of the system is re¬ 
quired to describe precisely the characteristics of the GW 
chirps, in particular the phase evolution. This may not 
be always possible as described in the following examples. 

The GW emitted by a coalescing binary of compact ob¬ 
jects can be divided into three phases (inspiral, merger 
and ringing). Although the GW can be obtained accu¬ 
rately in the inspiral phase when the bodies are well sep¬ 
arated [2 (using post-Newtonian expansions) and in the 
ringing phase after they have merged Q (using perturba¬ 
tive methods), the in-between merger phase still defeats 
both the numerical and analytical efforts 0 for modelling 
its highly non-linear regime. For large mass binaries, the 
merger phase contributes to a dominant fraction of the 
signal-to-noise ratio (SNR) Q. In this case, the search 
method has to accommodate the significant lack of signal 
information. 

Kerr black holes accreting matter from a surrounding 
magnetized torus areputative sources of the long gamma- 
ray bursts (GRBs) [g. It is claimed that, the black hole 
spin energy is radiated away through GWs along with the 
GRB. The precise shape of the emitted waveform would 
need accurate hydro-dynamical numerical simulations. 

A third example is the GW emitted in the form of 
the quasi-normal modes 0 by a newly born hot neutron 
star (during the cooling phase which follows the core col¬ 
lapse). Here, the characteristics of the GW depend on 
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the equation of state of the proto-neutron star and var¬ 
ious physical processes (like neutrino diffusion, thermal- 
ization and cooling) which are currently not known with 
accuracy. 

All these three examples are expected to emit GW as 
an unmodelled chirp, the phase information being not 
(perfectly) known. Its typical duration in the detector 
bandwidth is of the order of a few seconds. 

While matched filtering is a well-known and efficient 
detection technique when a precise waveform model is 
available, the lack of waveform information prevents us 
from using the same approach. It is thus natural to ad¬ 
vocate for exploratory searches (based on partial infor¬ 
mation or “good sense” models) as opposed to targeted 
ones (relying on a precise model). 

Various strategies Hi HI HI HI HI HI have been 
designed following this viewpoint, for the detection of 
transients of short duration (tenth to hundredth of mil¬ 
liseconds) or GW burst. Such transients are typically 
from supernovae core collapses. The notion of a varying 
frequency is not adequate for such a small number of cy¬ 
cles. It is thus not meaningful to describe such transients 
as chirps. Their detection is a different issue than the one 
considered here. 

Here, we are interested in exploratory search specifi¬ 
cally for unmodelled chirps. In the past, this question has 
already been investigated yielding a detection method, 
the Signal Track Search ^ (STS). The STS relies on 
the observation that, in a time-frequency (TF) represen¬ 
tation, a chirp appears as a filiform pattern and this dis¬ 
criminatory signature can be searched for. A satisfactory 
implementation of this phenomenological argument calls 
for a proper TF representation (TFR) and pattern search 
algorithm. The STS results from “ad-hoc” choices for the 
above mentioned points. 

In this paper, we propose a new method for the de¬ 
tection of unmodelled chirps. It is based on the same 
general principles (pattern search in a TFR) as the STS. 
Its originality resides in the clear link we establish be¬ 
tween the method (i.e., the choices of TFR and pattern 
search) and an optimality criterion. 

The paper is organized as follows. We state the detec¬ 
tion problem in Sec. im We introduce the general chirp 
model referred to as smooth chirp and we assume that 
most physically realistic GW chirps can be described by 
this model. The phase of a smooth chirp is an arbitrary 
continuous and differentiable function with bounded first 
and second derivatives. In Sec. nm we derive the opti¬ 
mal statistic for detecting a given smooth chirp in noise, 
which is usually referred to as quadrature matched fil¬ 
tering. The idea is then to apply this statistic for any 
smooth chirp, and select the maximum which is asso¬ 
ciated to the individual that best fits the data. This 
maximization has to be done numerically. To do so, the 
set of smooth chirp being a continuous set, has to be dis¬ 
cretized. In Sec. lEl we show that grids of templates can 
be constructed for smooth chirps using chains of small 
chirps, we call chirplet chains (GG). We further prove 


that the grid is tight i.e., any smooth chirp can be closely 
approximated by a chirplet chain. The maximization of 
the statistic over the set of smooth chirp can be reliably 
replaced by a maximization of the set of GGs. However, 
the number of GGs being very large, the computation of 
the quadrature matched filter for all GGs is not tractable. 
In Sec. m we propose a feasible (TF based) procedure 
for finding the best GG. We show that the quadrature 
matched filter can be reformulated approximately as a 
path integral computed in the TF representation given 
by the discrete Wigner-Ville (WV) distribution. As a 
result, the maximization of the statistic over the GGs 
amounts to obtaining the TF path of largest integral. We 
demonstrate that this kind of problem can be solved efh- 
ciently with dynamic programming. We detail our path 
search algorithm and we evaluate its computational cost. 
Finally, we compare the resulting algorithm with other 
methods in Sec. EH Receiver operating characteristics 
obtained in several realistic situations demonstrate the 
superiority of the proposed approach. 

II. SMOOTH CHIRPS IN GAUSSIAN NOISE 

We introduce a general chirp model which we refer to 
as smooth chirp, 

s{t) = Acos{(j){t)-\-tpo) for to < t < to + T’) (1) 

and s(t) = 0 outside this interval. 

A smooth chirp is characterized by the amplitude A, 
the initial phase (po and a smooth phase evolution (j)(t) 
(without loss of generality, we assume (j){t) = 0 at the ar¬ 
rival time t = to). We define the term smooth as follows. 
A phase (j){t) is smooth if this function and its first three 
derivatives are continuous and we have 



< F 


dt 


dt'^ 


for all t and where f{t) = {2 tt) ^dfi/dt is the instanta¬ 
neous frequency. The chirping rate limits F and F are 
chosen based on the allowed upper bounds obtained from 
general astrophysical arguments on the GW source of in¬ 
terest. The chirp model thus includes four parameters 
p = {A, ipQ,to, (j){-)} which are not known a priori and 
need to be determined from the data. 

Let the signal be correctly sampled at the Nyquist rate 
fs = 1/ts and let assume that we acquire the data Xk by 
blocks of N samples. The signal is denoted by Sk = s{kts) 
for fc = 0,..., A^ — 1 with the duration T = tgN. The 
noise is assumed to be additive white and Gaussian 
with zero mean and unit variance. Since the noise of 
GW detectors is colored, this noise model implies that 
a whitening procedure has been already applied to the 
data. (Therefore, the signal Sk in Eq. ((3J) is a “whitened” 
version of the actual GW signal). 

In this initial work, we restrict the smooth chirp model 
to have a constant envelope, although GW chirps are gen¬ 
erally amplitude modulated. The constant envelope thus 
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limits the descriptive power of the model. However, we 
argue that the model is still reasonable for many cases^ 
and that the phase information plays a major role for 
detection of chirps. We leave the problem of detecting 
amplitude modulated chirps for subsequent work. 


III. OPTIMAL DETECTION OF A SMOOTH 
CHIRP 

For each block of N data samples, the signal detection 
problem is to decide which statistical hypothesis suits 
best to the data among the following two : 

(iJo) Xk = rik noise only (3a) 

{Hi) Xk = Sk + nk signal+noise (3b) 

In practice, this requires thresholding a functional of 
the data, commonly referred to as statistic. If the statis¬ 
tic crosses the threshold. Hi is chosen as opposed to Hq 
and vis a versa. 

Due to the presence of random noise, this decision may 
not be always the right one. There are two types of er¬ 
rors associated to this: false alarms (decide Hi while 
Hq is present) and false dismissals (the opposite). The 
probabilities of occurrence of these two errors fully quan¬ 
tify the performance of a given statistic. This informa¬ 
tion can be used to rank the large number of possible 
statistics and to identify the best one. This is the ap¬ 
proach followed by the Neyman-Pearson (NP) criterion 
|l5| : the NP-optimal statistic minimizes one error prob¬ 
ability, while keeping the other fixed to a given value. 
To be precise, in the present case, it minimizes the false 
dismissal probability for a fixed false alarm probability. 

For simple detection problems, it can be shown 
that the likelihood ratio (LR) defined by A = 
F{{xk}\Hi)/F{{xk}\Ho) is NP-optimal ^3. For smooth 
chirp detection problem described in Eq. the LR can 
be easily obtained if we assume that the chirp parameters 
p are known in advance. When the parameters are not 
known a priori (which is the situation here), the ideal 
would be to have a statistic which is NP-optimal for all 
values of the parameters. This statistic is usually referred 
to as uniformly most powerful. However, it is not guaran¬ 
teed that such statistic always exists, and even if it does, 
it is generally difficult to obtain. 

A sensible solution consists in getting some kind of es¬ 
timates for the unknown parameters and then use the LR 
assuming that the estimated value is the actual value. If 
we use maximum likelihood (ML) estimators of the un¬ 
known parameters, the resulting statistic is referred to 


^ It is important to stress here that the model applies to the 
“whitened” chirp. For inspiralling binary chirps crossing the 
entire detector bandwidth, the envelope of “whitened” chirp is 
flatter than the original GW signal. For the other cases, this fact 
depends on the location of the chirp within the detector band. 


as generalized likelihood ratio test (GLRT) M (or max¬ 
imum likelihood test in the statistical community). 

The GLRT can be shown to be uniformly most pow¬ 
erful in certain cases M For our problem, up to our 
knowledge, this is an open question. Strictly speaking, 
it is thus not correct to qualify the GLRT as “optimal” 
(as is often done in the literature on GW data analysis). 
Nevertheless, we continue this misuse of language since 
the GLRT has proven to perform reasonably well and no 
better alternative appears to be available. 

In the following subsections, we give the derivation of 
the GLRT statistic. We proceed with the maximization 
of likelihood ratio with respect to the parameters. Fol¬ 
lowing we note that out of the four parameters. A, 
ipo and to are extrinsic parameters (known as kinemati- 
cal or dynamical parameters) whereas 4>{-) is an intrinsic 
parameter (which determines the shape of the chirp wave¬ 
form). On the basis of this distinction, the maximization 
over the extrinsic parameters can be treated in a sim¬ 
ple manner whereas the computation of the ML estimate 
of the intrinsic parameter requires a more sophisticated 
numerical treatment. 


A. Maximize the likelihood ratio: A and ipo 


In this subsection, we maximize the LR with respect to 
A and (po. In case of Gaussian noise, it is more convenient 
to use log-likelihood ratio (LLR) which is expressed by 


N-l ^ Af-1 

A(a;; p) = In A = ^ XkSk “ 2 ^ 


(4) 


We introduce Sk = cos((()fc -I- ipo) (such that Sk = Ask) 
with the norm Af = ^1- 

The maximization of the LLR A(a;; p) over A is 
straightforward and gives the expression of the ML es¬ 
timate of the amplitude, namely A = '^^=0 
Inserting this expression into Eq. 0 , we obtain 


A{x;{A,ipo,to,(l)i-)}) 


1 

w 




(5) 


The analytical maximization of the LLR over ipo de¬ 
serves a little more attention. The same calculation has 
been performed for the detection of chirps from inspi¬ 
ralling binaries min but it is based on the assumption 
that Af is independent of ipo which is not valid in the con¬ 
text of arbitrary chirps. In Appendix El we detail this 
calculation and discuss the validity of this assumption. 

We express the resulting statistic £{x]to,(j)) = 

A(a;; {A, (^0: ^0, <('(•)}) using the following notations for 
the cross-correlation of the data with the two quadrature 
waveforms. 


W-l N-l 

Xc= Xk COS ct)k Xs= Xk sin (l)k , (6) 

fe =0 / c =0 
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and for the norms and cross-products of cos (j)k and sin (j)k , 


N-l 


N-l 


E 

cos^ (j)k 

n* = ^ sin^ (j)k 

(7a) 



k^O 



N-l 




-E 

COS (/)/c sin (/)/e . 

(7b) 


/c=0 




We distinguish two cases. In the degenerate case where 
the two quadrature waveforms are linearly dependent ((j)k 
is a constant), O = ncTis — vanishes and we have 

e{x;to,(l>) = {xl+xl)/{2N). (8) 

Otherwise O > 0, the optimal statistic is 

£{x; to, (j>) = {risxl - 2na;XcXs -f 71^x1)1(20), (9) 

and is commonly referred to as quadrature matched fil¬ 
tering, (see Appendix IXll. 

B. Maximize the likelihood ratio: (j) and to 

The statistic i results from a quadratic combining of 
the cross-correlations defined in Eq. ©■ It can be seen 
as a “generalized dot-product” and can be related to a 
“distance” measuring the discrepancy between the data 
and template waveforms (or, in short, templates) defined 
by the phase (see Eq. (IA8II 1. Maximizing £ over cj) is 
equivalent to minimizing this distance. 

The expression in Eq. is for a given known phase fi. 
If the phase is unknown but belongs to the set of smooth 
chirps, then we need to minimize the distance within this 
feasible set. In other words, we need to find that smooth 
chirp which best fits the data i.e., find 

£inaAx;to) = max {£(x; to, ()>)}■ (10) 

all smooth chirps 

This maximization is difficult to tackle analytically and 
has to be done numerically. The set of smooth chirps 
is a continuous set and hence not easy to manipulate 
numerically without discretizing it. For this purpose, we 
introduce chirplet chains, which we discuss in the next 
section. 

As described earlier, we process the data stream block- 
wise. We compute the statistic independently for each 
block. The maximization over tg is obtained by compar¬ 
ing ^max for neighbouring blocks and selecting the maxi¬ 
mum. The ML estimate of to is then given by the starting 
time of the corresponding block. The period separating 
two successive starting times thus defines the resolution 
of the estimate. If required, this resolution can be im¬ 
proved by increasing the overlap between two neighbour¬ 
ing blocks. 

We now concentrate on the maximization of £(x', (j)) 
over 0 in a given block. In the following, we remove to 
from the arguments of £ to keep the notations simple. 


IV. CHIRPLET CHAINS: A TIGHT TEMPLATE 
GRID FOR SMOOTH CHIRPS 

In this section, we show that chirplet chains (CCs) can 
be used to construct template grids for smooth chirps. 
CCs are based on the simple geometrical observation: 
broken lines give good approximations of smooth curves. 
CCs are signals whose (instantaneous) frequency is a bro¬ 
ken line. We verify that they are good approximation of 
the frequency curve of an arbitrary smooth chirp. We ob¬ 
tain the conditions ensuring that, for any smooth chirp, 
there always exists a sufficiently close CC. If these condi¬ 
tions are satisfied, the set of the CCs forms a tight tem¬ 
plate grid which can be used to search for an unknown 
smooth chirp. Finally, we examine the implementation 
of such grid for the toy (but realistic) model given by the 
inspiralling binary chirp. 


A. Chirplet chains: piecewise linear frequency 

All smooth chirps in Eq. o are supported in the TF 
domain V, a rectangle of width T and of height equal to 
the Nyquist bandwidth /s/2, as illustrated in Fig. ^ Let 
{(tj = jSt, fm = mSf) ; j = 0...Nt, m = 0... Nf} be a 
regular TF grid led on T> by splitting the time axis into 
Nt intervals of size 6t = T/Nt, and the frequency axis 
into Nf bins of size 5f = fs/{2Nf). 

In the following, the subscripts j and m designate the 
index of the time interval and the frequency bin of the 
grid, respectively. The index k G {0,..., A — 1} denotes 
the time index of a sample. 

A chirplet is a short piece of signal whose frequency 
varies linearly between two successive nodes of the grid. 
In the time interval j, we denote the time and frequency 
coordinates of the chirplet extreme points by (j, rrij) and 
(j -\- I,TOj+i). In the TF plane, it is thus represented by 
a line joining the grid nodes (tj,fmj) and (tj+i,fmj+i) 
(see Fig. Concretely, this means that the phase fik = 
4>(tsk) of a chirplet is a quadratic function of time, as 
follows, for tj < kts < tj+i 

— ^j£j,k ffi ^j^j,k T — 

where aj = ir (/mj+i - f'mj)/dt, bj = 27r and tj^k = 

tgk — t j . 

We build the chirplet chain (CC) by enforcing chaining 
rules. The frequency and phase of this chain are contin¬ 
uous. Clearly, the continuity of the frequency is ensured 
by construction, while the phase continuity requires that 

^i-l = + fnij-i) + dj-2 , (12) 

for j > 1, and fixing 6*_i = 0. The slope of the 
chirplet frequency as well as the difference between the 
slopes of the frequencies of two consecutive chirplets are 
bounded absolutely. These bounds are given by the 
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FIG. 1: Chirplet chains — The TF domain of interest 
T> is tiled by a regular TF grid of Nt time intervals and Nf 
frequency bins. A chirplet is a short piece of signal whose 
frequency varies linearly between two successive nodes of the 
grid. It is thus represented by a line joining the grid nodes. 
The slope of the chirplet frequency is limited (triangular re¬ 
gion in light gray, here N'^ = 1). We chain the chirplets, 
imposing the continuity of the chain and limiting the differ¬ 
ence between the slopes of two consecutive chirplets (triangu¬ 
lar region with dark gray stripes, here N” = 1). Admissible 
chirplets in time interval j belong to the intersection of these 
two regions associated to the regularity constraints. Clearly, a 
chirplet chain is represented by a broken line in the TF plane. 


two parameters N'^ and TV" respectively such that (z) 
\mj+i — rrijl < and {ii) \mj+i — 2mj + < TV". 

Clearly, a CC is represented by a broken line in T>. 
The two parameters and N” control the regularity 
of this line. Consistently, we will refer to (z) and (ii) as 
regularity constraints. 

The instantaneous frequency of a smooth chirp is asso¬ 
ciated to a smooth curve in T>. In the same manner that 
broken lines are good approximations of smooth curves, 
CCs are good approximations of smooth chirps. Since 
CCs form a finite discrete set, they sample ^ the set of 
smooth chirps. In other words, they form a template grid 
of this set. 

It is important to know whether this template grid is 
sufhciently tight i.e., whether for any smooth chirp, there 
always exists a sufficiently close CC. The template grid 
tightness is controlled by the choice of the four parame¬ 
ters defining the set of CCs, namely the TF grid param¬ 
eters Nt, Nf and the regularity parameters N^ and TV". 
The first and preliminary step to address the tightness 
question is to define a distance measuring the “similar¬ 
ity” (or ambiguity) between two different chirps. 


^ Strictly speaking, CCs don’t sample the set of smooth chirps 
since they don’t belong to this set (the second derivative of their 
frequency is not defined at the boundaries of the grid time inter¬ 
vals and it is thus not bounded). 


B. Distance in the set of smooth chirps 

We follow the approach suggested in 0 and assume 
that we “receive” a chirp whose phase (p is different than 
the template phase p*. We set Xk=Sk = Acos{(j)k + 9), 
and consider 


c{p,p*) 


£{s;(p) 


(13) 


Clearly, £ measures the reduction factor of the “detec¬ 
tion peak” due to the mismatch between the chirp present 
in the data and the chosen template. It is a relative mea¬ 
surement done with respect to the ideal case where the 
template matches exactly the considered chirp. In this 
sense, it can be interpreted as a SNR loss. 

Since £ > 0 and equals 0 when (p = (p*, it can be 
interpreted as a distance between the chirps. Note that £ 
does not depend upon A. It depends only on the phases 
p and p*, but this dependency is difficult to perceive 
intuitively from its definition in Eq. m- 

An approximated but much simpler expression can be 
obtained for nearby chirp and template retaining the 
leading terms of a Taylor expansion for small Aj- = 
Pk ~ Pk- The approximation is detailed in Appendix H 
and leads to the following expression 

1 ^-1 

(14) 

k^O 

with A = 1/NY,^~q Afe. 

Interestingly, we recognize in this expression the em¬ 
pirical estimate of the variance of the phase difference 
Afc. With this definition of the distance, two chirps are 
“identical” (their distance measured by £ is zero) if and 
only if they have the same phase evolution up to an ad¬ 
ditive offset. 


C. Is the CC grid tight ? 

In this section, we address the grid tightness problem 
and find the regularity and TF grid parameters which 
yield a tight template grid of CCs. We proceed as fol¬ 
lows: we first consider an arbitrary smooth chirp of phase 
p. Then we construct a CC “geometrically” close to this 
chirp. We check if this CC is admissible i.e., if it satisfies 
the regularity constraints. This imposes two conditions 
on the parameters. Finally, we check whether it is ef¬ 
fectively close to the chirp (as measured by the metric). 
This yields the loss due to the approximation of the chirp 
by a CC in the worst case. For tight CC grid, this loss 
(homogeneous to a SNR loss) has to be small which im¬ 
poses one more condition on the parameters. 
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arbitrary chirp 
"geometricaliy" dose CC 


FIG. 2: (Color online) Smooth chirp (solid cnrve) and its 
“geometrically” close CC (dotted broken line) 


We rewrite this condition in the following form 


N' Nf 

N' > 4—- - + 1 

Nt2N^ ^ 


(18) 


where N' = FT^ is an adimensional quantity which 
depends only on the fundamental characteristics of the 
smooth chirp model. 

b. 2"^^ order regularity — We consider two succes¬ 
sive chirplets in intervals j — 1 and j. Using a similar 
method, the difference of their slopes can be bounded by 


1. “Geometrically” close CC 

Let (j){t) and f{t) be the phase and frequency of an 
arbitrary smooth chirp. The frequency evolution appears 
as a smooth curve in the TF plane. 

We construct a CC “geometrically” close to this chirp 
as follows: for each j = 0,... ,iVt, we choose the node 
{tj, /*) of the j-th column of the TF grid defined in 
Sec. 11 V Al which is nearest to the point (tj,f{tj)). We 
draw the broken line by joining these nodes (see Fig. [2J|. 
The associated CC is the “geometrically” close CC to the 
chirp under consideration and we denote its phase (j)*. 


2. Admissibility of the “geometrically” close CC 

For a given chirping rate limits F and U, the “geomet¬ 
rically” close CC may or may not satisfy the regularity 
constraints. This depends on the regularity and TF grid 
parameters. Below, we investigate this question. 

a. 1^* order regularity — Let us consider the 
chirplet of the interval j, we have 

\fUi - /;i < 1/^1 - fih+i)\ + \fitj+i) - 

+ (15) 

Using the mean value theorem^ (see e.g., M), we get 
l/(0~/(s)l — from which we deduce a bound on 

\f itj+i)-f{tj)\- By construction, we have \ f*-f{tj)\ < 
Sf/2 and this leads to 


i/;vi - 2 /; + f;-i\<- 2f{t,) + /(t,_i)i + 25f. 

(19) 

Two consecutive applications of the mean value the¬ 
orem to a function /(•) which satisfies Eq. |2Jl for all 
s G [r, t] with 0 < r < t < T yield the following result 

\f{t) - 2/(s) -f /(r)| <F{{t- sf + {r- s)^) /2 

+ F\t-2s + r\. (20) 


Using r = s = tj and t = tj+i, we get 

\f*+,-2f; + f;_^\<F5l + 25f. (21) 

Therefore, the geometrically close CC satisfies the reg¬ 
ularity constraint {ii) mentioned in Sec. CVUif 


N'; > FS'^/Sf + 2 . 


We rewrite the above condition as 



/N" 


2 


Nf 

_ L 2 

2N ^ ’ 


( 22 ) 


(23) 


where N" = V3FT^ is an adimensional quantity which 
depends only on the fundamental characteristics of the 
smooth chirp model, and thus from the astrophysical in¬ 
put. It is related to the maximum overall curvature of the 
chirp frequency or more precisely to the largest number 
FT^ of Fourier bins that the chirp frequency can sweep, 
the linear trend being removed. 


\f;+,-f;\<Sf + F6t. 


(16) 


3. Is the “geometrically” close CC effectively close? 


Thus, the geometrically close CC satisfies the regular¬ 
ity constraint (i) mentioned in Sec. IIV Al if 


K > F^t/Sf + 1. (17) 


® Let the function g(-) be continuous in the open interval (a, b) 
and differentiable in the closed interval [a, b]. The mean value 
theorem states that there exists c in (a, b) such that gib) — g{a) = 
g(c)(6 — a) where g(-) is the derivative of g(-). Consequently, we 
have \g{b) - g{a)\ < G\b - a| with G = sup„,g(,, (g(a:)). 


We obtain a worst case estimate on the distance be¬ 
tween the smooth chirp and its “geometrically” close CC 
by bounding the variations of their phase difference. We 
begin with bounding the frequency discrepancy. 

The starting point is the following lemma (inspired 
from 1^, p. 23) obtained from the application of the 
mean value theorem and some algebraic manipulations. 
If /(•) satisfies Eq. |(2Jl for all s G [r, t) then 


f{s) 


f{r) + 7 —- f{r)) 

t — r 


< F{t 


rf. (24) 
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However, this upper-bound can be slightly improved: 
the term (t — r)^ over-estimates the more precise bound 
g{s) = min{(t -s)(t-{s + r)/2), (s - r)((s -h t)/2 - r)}. 
(Note that g{t) = 0 and g{r) — 0 as expected.) In the 
worst case, we have s = {t+r)/2 and g{s) = 3/8(t—r)^ as 
opposed to (t —r)^. We include this gain in the following. 

We apply the lemma in Eq. (I24II with t = tj+i and 
r = tj. The term inside the square brackets in the left 
hand side of Eq. is equal to the frequency at time s 
of the chirplet obtained by joining to f{tj) (see 

the dot-dash line in Fig. El. We denote this frequency 
/(s), and then obtain 

\f{s)-m\<3F6^j8. (25) 

^Since \ f{s)-f*{s)\ < \f{s)-f{s)\ + \f{s)-f*{s)\ and 
|/(s) - /*(s)| < Sf/2, we have 

\f{s) - ns)\ < 3FSf/8 + Sf/2 = Af. (26) 

By definition (j){t) — 4>{r) = 27r f{s) ds. Integrating 
both sides of the above inequality between two successive 
points r = ts{k — 1) and t = tsk for A: e {I,..., iV — 1}, 
we get 

|Afc-Afc-i| < 27rA/t«, (27) 

which constraints the variations of the phase difference 

Ak = r-4‘i- 

We prove in Appendix o that the approximated dis¬ 
tance C{(j),(j)*) as shown in Eq. 11411 is maximum under 
this constraint, when A^ = ±.27: A ft gk. In this case, the 
maximum is £{(/), tp*) = (7rAyT)^/3 = /j,'. 

We can finally state the following tight template grid 
theorem: for all smooth chirps of phase (p, there exists a 
CC of phase (p* such that 

C{<P,cP*)<ti', (28) 

where g' = {3FS^ / 4 + 6 f)'^ /12 is the maximum (i.e., 

in the worst case) energy SNR loss due to the mismatch 
between the smooth chirp and the chosen template given 
by a close CC. The corresponding maximum amplitude 
SNR loss is ^ = 1 — ~ ^ j2 for small g!. Note 

that the amplitude SNR loss is linked to the minimal 
match MM defined in by the relation MM = 1 — g. 
We express the maximum SNR loss g in terms of the CC 
parameters as 


TT^ r 1 1 /2A\ 


(29) 


In principle, this loss can be made arbitrarily small by 
choosing Nt and Nf adequately. Therefore, the grid of 
CC can sample the set of smooth chirps tightly. 

It is evident that two types of losses contribute to g. 
The first one is related to the geometrical error due to 
the fact that the model is a broken line: within the 


time intervals of the TF grid, the model is a straight 
line which cannot perfectly follow the curvature of the 
smooth chirp frequency. The finer the grid along the 
time axis, the smaller the time interval, the better the 
line hts the smooth chirp frequency, thus reducing this 
error. The other is related to the quantization error as we 
require the node of best broken line to belong to the TF 
grid: there is a difference between the best broken line we 
can possibly draw and the closest (quantized) one with 
vertices belonging to the grid. The hner the grid along 
the frequency axis, the closer the quantized line from the 
original, thus reducing this error. The maximum SNR 
loss is the function of these two independent parameters. 

When Nt = N" and Nf = 2N , the maximum SNR loss 
is of order ^ 7r^/96 « 10% and the two types of errors 
contribute equally. The same maximum SNR loss can be 
achieved with other choices for Nt and Nf. In the next 
section, we propose a criterion to solve this indetermina¬ 
tion. 


D. Smallest tight CC grid 

As elaborated in the previous section, the tight tem¬ 
plate grid theorem gives the condition on the TF grid pa¬ 
rameters Nt and Nf which ensures that the template set 
(the CCs) covers all the feasible set (the smooth chirps) 
with a given accuracy specihed by the maximum SNR 
loss. The same accuracy can be achieved with several 
pairs of parameters, leading to a parameter indetermina¬ 
tion. 

For a small maximum SNR loss, the maximization of 
the LLR in Eq. oni performed over the set of smooth 
chirps can be safely replaced by a maximization over the 
set of CCs, i.e., 


■^max (cc) ~ max {£(a;; ()))}. (30) 

all CCs 

The statistic ^max in Eq- results from the CC 
which maximizes the statistic or in other words, from 
the waveform of the template set which best fits the 
data. Generally speaking, when the data is noise only, 
the larger the number of (reasonably different) waveforms 
in the template set, the larger the risk that one of the 
waveforms fits the noise and consequently, the larger the 
false alarm rate. 

The TF grid parameters Nt and Nf influence very dif¬ 
ferently the number of CCs. The above argument sug¬ 
gests to select the parameters which minimize the num¬ 
bers of CCs, for a given specified maximum SNR loss. 
We refer to the smallest tight CC grid as the set of CCs 
which results from this constrained optimization. 

Let us first estimate the number of CCs. According 
to the regularity conditions, each of the number Nc ^ 
Nf{2NP + 1) of possible chirplets in a given time interval 
can be chained to (at most) 2N" -|- 1 chirplets in the 
next time interval. Counting CCs is then a combinatorial 
problem. We have Nc chirplets in the first time interval. 





and 2iV" + 1 possible choices for the Nt — 1 successive 
time intervals. Neglecting what happens at the lower and 
upper boundaries of the frequency axis (i.e., near DC and 
Nyquist), we obtain an upper-bound on (the logarithm 
of) the number Ncc of CCs as 


We consider the Newtonian approximation of the chirp 
whose frequency evolution is given by [T^ 

/ f _ t \ 

fit) = /o f 1- jAj for t<to+T, (35) 


In Wc < H^KNf) + (Nt - 1) H‘2K' + !)• (31) 

In practice, we have Nt ^ 1. The second term largely 
dominates the right hand side and the first term can be 
neglected. We thus have IniVcc ^ Nt \n{2N” + 1). 

At this point, it is convenient to introduce u = N"/Nt 
and V = 2N/Nf and express the smallest tight CC grid 
problem with these variables. From the regularity con¬ 
straints, we have N” = 4m^/(3u) -I- 2. We want to mini¬ 
mize the number of CCs 

InNcc oc g{u,v) =-\n (^—+ 5 ] , (32) 

U \6 V J 


subject to a given maximum SNR loss i.e., vf + v = C = 

Sy/^/lT. 

Combining the derivatives of the objective dg = 
dugdu -|- dygdv and of the constraint dv = —2udu, we 
obtain the equation giving the admissible point where 
the derivative dg/du vanishes, viz. 


y - 5 


1-^=0 
4y 4 


(33) 


where we defined y = 8u^/(3u) -I- 5. This equation 
can be solved numerically and gives y ~ 8.95. Let 
a = vf/v = 3{y — 5)/8 be the ratio between the two 
errors contributing to y. We obtain the smallest tight 
template grid when this ratio is a ~ 1.48. For a required 
/i, we get the parameters of the resultin g grid as foll ows. 
Using the constraint, we have u = Caj (1 -I- a) and 
V = C/(1 -\- a), from which we obtain the parameters, 


Nt = 0.52 N" , Nf = 0.78 ( 34 ) 


Interestingly, this also implies that N!/ = 4a/3-|-2Ri4 
is a constant (i.e., does not depend on y). The last 
parameter N/ is directly determined by substituting 
Eqs. (IS3|) in (IT^ . 

The parameters of the smallest tight template grid may 
not be always suitable in practise (see the later discussion 
on the implementation and numerical contingencies in 
Sec. IWl but they give interesting indications. 

At this point, it is useful to see with an example if the 
proposed model and template grid sound tractable in a 
realistic case. 


where to denotes the arrival time. In practice, the arrival 
time corresponds to the time at which the chirp enters 
the detector’s bandwidth i.e., when its frequency reaches 
the low frequency (seismic) cut-off (denoted /o) of the 
interferometric detectors. The T defines the chirp dura¬ 
tion, i.e. the time taken by the chirp from the arrival 
time till the binary coalescence. 

The chirp duration can thus be estimated by 


T 



M 

50Mq 



(36) 


where M is the total mass (objects of equal masses). 

In this calculation, we assume the seismic cut-off fre¬ 
quency"* of 20 Hz. 

We fix F and F to the corresponding values of the first 
and second derivatives of the chirp frequency, pertaining 
to the last stable circular orbit (LSCO^) viz., 

/lsco~88.4Hz(^^^ , (37) 

""-""'“’’/“'(ail;) 


We note that T, /lsco, F and F decrease with an in¬ 
creasing mass. When M increases, the chirp is thus 
shorter, less steep and curved, and it reaches only the 
lower part of the frequency band. From the above equa¬ 
tions, we deduce that 

M 

50Mq ) 
(40) 

The sampling frequency fs is fixed by the width of 
the observational band of the GW detector, namely fs = 
2048 Hz. We thus have 


N' ~ 2.2 X 10^ 


/ M 

V5OM0 


-16/3 


N" 


698 


/ M 

N = fsT 2662 - . (41) 

^ \50MqJ ^ ’ 

Following Sec. MB and fixing y = 10 %, the smallest 
tight CC grid has the following parameters for the TF 


E. Toy model and CC parameters 


We use the inspiralling binary chirps as a toy model 
to check whether the various parameters have reasonable 
order of magnitudes in this physically realistic situation. 


This is the seismic cut-off frequency targetted by the detector 
Virgo. 

® For non-rotating stars, the LSCO is when the objects are at the 
distance r = etGMjf?. 
















9 


grid 


B. Near optimal search 


Nt ~ 645 


Nf ~ 6566 


M 

50MfT 


-4 




M 


\50Mq 

and for the regularity, we have 


-5/3 


n: -17 


M 


50Mq 


-4/3 


NL' - 4. 


(42) 

(43) 

(44) 


The orders of magnitude for the various parameters ap¬ 
pear to be reasonable. Since these parameters don’t in¬ 
crease with M, the template grid defined with the above 
values remains acceptable and tight for higher masses 
M > 50Mq. 


The maximization of £{x] (p) in Eq. (I30II is a combinato¬ 
rial maximization problem. The existence of an efficient 
solving algorithm for such problem is related to the struc¬ 
tural properties of the “objective” function to be maxi¬ 
mized, that is, £ in the present case. In this Section, we 
show that £ can be reasonably approximated by a path 
integral computed over a time-frequency representation 
(TFR) of the data. The structure of the approximated 
statistic allows us to perform its maximization efficiently 
with dynamic programming. The approximation goes 
through two stages with an intermediate step for the re¬ 
formulation of the statistic in the TF plane. 


1. Approximation 1: for a CC, cosine and sine are almost 
orthogonal 


V. FIND THE BEST CHIRPLET CHAIN 


As shown in Eq. HASH , the statistic £ can be expressed 
as 


In Sec. nvi we have shown that, the SNR loss due to 
the use of a CC instead of the ideal template can be made 
small with an appropriate choice of parameters i.e., by 
making the CC grid tight. In other words, the problem 
of detecting a smooth chirp is equivalent to the one of 
detecting a CC as stated by Eq. (EOl)- The maximization 
over the set of CCs - involved in the latter case - has the 
great advantage that it can be resolved numerically. 


A. The exhaustive search is not feasible 

Since CCs are in finite number, an obvious maximiza¬ 
tion procedure is to try them all and select the one which 
gives the maximum. To understand whether this solution 
is tractable, we need to know how many CCs are there. 
We consider that the search parameters At, Ay, A' and 
N" are known and can be obtained from the physical and 
grid tightness requirements as discussed earlier. 

We already presented an estimate of the number of 
CCs in Eq. (I3III and saw that it grows exponentially 
with the number of time intervals of the TF grid. This 
estimate computed for the toy model example presented 
in the previous section gives log^Q Ncc ~ 1400. Clearly, 
this number is too large for an exhaustive search (i.e., 
computing £ for all possible CCs) to be carried out in real 
time on existing computers. Generally speaking, since 
the number of CCs increases exponentially with A*, the 
cost of an exhaustive search scales exponentially with At 
and thus with the problem size A. 

In the next section, we propose an algorithm which 
gives a good estimate for the optimal CC instead of the 
exact solution of the maximization problem described in 
Eq. 13011 . However, as opposed to the exhaustive search, 
the computational cost of this algorithm scales as a poly¬ 
nomial of the problem size A. 


(x; (p) = 


I 


/N-l \ /N-1 \ 

^ I T f ^ ) ^k^k 

\ k=0 ) \ k=0 / 


(45) 


where the templates Cfc and are the orthonormalized 
counterparts of the waveforms in quadrature cos (pk and 
sin (pk obtained from the Gram-Schmidt procedure as 
given below 


cos (pk 



Sk 


ric sin pk - rix cos pk 
VricO 


(46) 


Let {cos(/)fc} and {sint/ffc} be the vectors in associated 
to the quadrature waveforms. As it appears in the above 
expressions, these vectors are generally not orthonormal. 
Their deviation from orthonormality can be quantified 
with two parameters, defined by 


ric - ris 
ric + ns' 


‘Ifix 

ric + ns' 


(47) 


which are related to their vector lengths Uc, Us and their 
scalar product rix- The parameter 5 measures the relative 
difference in the vector lengths while e measures the angle 
between them. The vectors are orthonormal if and only 
if both 5 and e are zero. 

Intuitively, if the quadrature waveforms oscillate suffi¬ 
ciently, they should be close to orthonormality and 5 and 
e are expected to be small. This intuition is examined in 
details in Appendix^ in which we exploit the fact that 
(p is not arbitrary but it is the phase of a GG. We show 
that if (p is the phase of a GG whose node frequencies are 
in the bandwidth fi < fm^ < fs/2 — fi for all j with 



( 48 ) 


then |i5| < 77 and |e| < 77 . 
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In the following, we assume that this condition is sat¬ 
isfied. This imposes the CC frequency not to approach 
arbitrarily close to the DC nor to the Nyquist frequen¬ 
cies. Since the amplitude of the instrumental noise of 
GW interferometers diverges rapidly when going close to 
DC, it is not expected to detect GWs at low frequencies. 
Therefore the reduction of the bandwidth in the low fre¬ 
quency region should not be a problem as long as fi re¬ 
mains small. We will check later with examples that the 
reduction of the useful bandwidth is indeed sufficiently 
small. 

Using Eqs. TO and TO . we can write Ck and Sfe in 
terms of S and e as 

/ 2 

f 2 (1 + 5)sin0fc - ecos^fc 

= (.vanjj 

noting that Uc = N{1 + 6)/2, Ug = N{1 — 6)/2 and Ux = 
iVe/2. 

Inserting this expression in Eq. TO and taking the 
limit for small <5 and e, we find that £(x; </>) ^ £{x] 4>) = 
{x1-\-xl)/N, the reminder R{5, e) = i{x; (j))—£{x; (j)) being 
given by 


R{5, e) 


1 5{x1 — xD + 2eXcXs — ((5^ -I- e‘^){xl + x'^) 
N 1 - (52 + £ 2 ) 


(51) 

Considering that we have ^2 -|- e 2 < Appendix 

rnt and 


\xl - a; 2 | 


x^+x"^ ~ 


< 1 , 


1 2XcXs 


< 1 , 


xt: 


(52) 


the relative error can be bounded as 

I^(^U)I ^ 277-772 

< --^ « 277 , 


£{x;(l)) 


1 — 


(53) 


for small 77 . 

Provided a good choice of 77 (and checking the con¬ 
sequences on fi), this approximation error can be made 
small. We can safely replace i hy £ which we express as 
the following complex sum : 


the TF domain, provided the use of a unitary TER. One 
such TER is the discrete Wigner- Ville (WV) distribution 
defined in and given by 

kn 

wx{n,m)^ £ (55) 

k— — kn 


with kn = min{2n, 2iV — 1 — 27i}, pn^k = [tt -I- fc/2j 
and Qn^k = [n — k/2\ where [-J gives the integer part. 
The arguments of Wx are the time index n and the fre¬ 
quency index m which correspond in physical units, to 
the time = tgU and the frequency is fm = fs'm/{2N) 
foi 0 < m < N and fm = fs{N — m)/{2N) foi N + 1 < 
m < 2N — 1. Thus, the frequency axis gets sampled 
at twice the usual rate (as performed by the EFT). The 
WV distribution is associated with a particular sampling 
of the TF plane. As discussed later in Sec. IVTT^ this 
leads to some restrictions on the TF grid used for defin¬ 
ing CCs. 

Let {a;^} and {yk} be two time series. Moyal’s formula 
states that 


yk 


^ N-12N-1 

^ wx(n, 777 ) 777 ^( 7 ^, 777 ). (56) 

n=0 m—0 


Using Eqs. iTO and we rewrite £ as the inner- 
product of two TFRs namely, the WV of the data Wx 
and the template WV We which is the WV of complex 
template waveform Ck = expi(j)k, 

^ A7-12Ar-l 

X! Wx{n,m)we{n,m). (57) 

n—0 m—0 

Chirp signals are easily modelled and described in the 
TF plane. Qualitatively, we expect that the TER of a 
chirp signal have large values essentially in the vicinity 
of a curve corresponding to their instantaneous frequency 
and vanishes elsewhere. The template WV We being the 
TER of a chirp, it shares these characteristics. In the 
following section, we make use of this feature to simplify 
the statistic. 


3. Approximation 2: the WV of a CC is almost Dirac 


£{x] (j)) 


1 

V 


Xkeyip{i(j)k) 

k^O 


(54) 


2. Go to time-frequency: Moyal 

The expression of £ in Eq. TO computes the canon¬ 
ical hermitian scalar product between the data and a 
complex template waveform. While Parseval’s formula 
allows an equivalent formulation of this scalar product in 
the frequency domain, Moyal’s formula does the same in 


With continuous time and frequency variables, it is 
well-known that (H^, p. 130 and also 217) the WV of 
a linear chirp (i.e., a chirp whose frequency is a linear 
function of time) is a Dirac distribution along the TF 
line associated to the chirp frequency. 

We assume that this remains reasonably true for dis¬ 
crete time and frequency and when the chirp is non-linear 
(and in particular, when it is a CC). More precisely, we 
consider that we have 

We(n,m) K 2N 5{m — rrin), (58) 

where 777 „ = [ 2 r/„] where [•] denotes the nearest integer. 
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Here, /„ is the instantaneous frequency of the (possi¬ 
bly non-linear) chirp. Eq. (I58II dissembles two approxi¬ 
mations which we explain now. 

For discrete time and frequency, the discrete WV of a 
linear chirp can be calculated analytically |^. For the 
positive frequencies i.e., for 0 < m < N, the model in Eq. 
(P|l is an acceptable approximation of the exact result, 
as illustrated in Fig. 0 However, there is a significant 
difference in the negative frequencies i.e., for iV -|- 1 < 
m < 2N — 1. In this region, the discrete WV exhibits 
aliasing terms (clearly seen in the left panel of Fig. |31) 
which are closely related to the unitarity property of the 
WV. In 1^ , the aliasing terms are shown to be oscillating 
terms (switching signs) with smaller amplitude than the 
preponderant terms modelled by Eq. We can then 

expect their contribution to the summation in Eq. 0 
to be negligible. 

It is well known |22|| that interference terms appear 
when computing (both continuous and discrete) WVs of 
non-linear chirps. They can be related to the quadratic 
nature of this distribution (see for a detailed analysis 
of the nature and geometry of these interferences). Inter¬ 
ference terms change sign rapidly (see Fig. 13 right panel) 
and can be neglected for the same argument invoked for 
aliasing terms. 

Inserting Eq. into Eq. we get the following 

approximation of I : 

1 

n—0 

We see that this statistic results from the integral of 
the WV of the data along the TF path determined by 
the CC frequency /„. In other words, this integral is the 
area under this TF path. We refer to this quantity as the 
yath length 

With this approximation, the maximization of the 
statistic in Eq. amounts to hnding the path giving 
the largest integral, or the longest path. Efficient meth¬ 
ods exist for longest path problems These methods 
exploit the structural properties of path length (or inte¬ 
gral) measurement, in particular, additivity. The length 
of this entire path can be measured by splitting the path 
and summing the length of its constituent parts. Thanks 
to this property, the maximization problem can be de¬ 
composed into a recursive series of small problems, each 
of them being solvable in polynomial time. This is the 
main principle of dynamic programming (DP), which we 
describe in the next Section. 

We remind the reader that, contrarily to the new 
statistic i, the exact statistic I is not additive. DP cannot 
be applied to maximize £. 


If we see the WV as a Lebesgue measure (although this is an 
misuse of language since the WV can take negative values), the 
integral in Eq. effectively defines a path length. 


linear non-linear 




FIG. 3: Discrete Wigner-Ville of two chirp signals — 

The signals are normalized to unit £2 norm and we show the 
contour at the level 1/8. left: when the chirp frequency is 
a linear function of time, its WV is almost Dirac along the 
corresponding TF line in the TF half plane associated to pos¬ 
itive frequencies. For the negative frequencies, the WV distri¬ 
bution exhibits aliasing terms (we highlight them by a back¬ 
ground in light gray) which we neglect in the simplified model 
in Eq. 1581 . right', when the chirp frequency is not linear 
(here, it is a parabolic chirp), interference terms appear (ev¬ 
idenced by a dark gray background). Their contribution are 
also disregarded in the simplified model. Note that the WV 
of the non-linear chirp chosen for this illustration do present 
aliasing terms (with a light gray background), but they have 
a smaller amplitude than in the linear case. 

4 . Maximization with dynamic programming 

DP is a classical method for solving combinatorial 
optimization problems. As explained in the previous sec¬ 
tion, the idea is to decompose the problem into smaller 
ones that can easily be solved. In our context, the natu¬ 
ral decomposition is given by the tiling of the time axis 
into chirplet intervals i.e., tj < tgU < tj^i or equivalently 
jh < n < (j -I- 1)6 — 1 where b = dt/tg is the number of 
samples in an interval. The overall path integral is equal 
to the sum of the integrals computed in each chirplet 
interval marked with an superscript index as follows 

Nt-l 

^ l'■^(a;;(/)) (60) 

j=o 

^ U+l)b-l 

with f {x; (j)) = — ^ Wa:{n,mi), (61) 

n=jb 

where = [2N fi] and the frequency // follows the line 
joining the grid points {tj.fm^) and {tj+i,fmj+i). We 
also denote with a subscript index j, the path integral 
up to interval j, viz. 

3 

^ [x',(j)). (62) 

j'=0 

DP relies on the principle of optimality. We elaborate 
this principle with the help of Fig. 0] We consider the 
chirplet in time interval j. In a chain passing through 
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this chirplet, the regularity constraints limit the choice 
of preceding chirplets in the time interval j — 1. We 
suppose that there are only three such chirplets; namely 
a, [3 and 7 . 

Now, consider the time interval j — 1. We assume that 
we know the chain passing through the chirplet 0 {z be¬ 
ing either a, jS or 7 ) and giving the largest path integral 
summed up to the interval j — 1. We denote this quantity 
by ^j-i- (In this discussion, the chirplet and its associ¬ 
ated CC are designated by the same label). 

We compute the path integral contribution in j-th in¬ 
terval for the considered chirplet, and add the result to 

^(z) '"(z) 

(-jli to obtain t- ’ for all the three paths z = a,(3 and 7 . 

We mark with (*) the optimal chain associated to the 
global maximum of £ (i.e., summing from interval 0 to 
Nt — 1) which we denote We further assume that 
this optimal chain (*) follows (a) up to interval j — 1 , 
continues following the considered chirplet in interval j 
and proceeds to the last interval j = Nt — 1 with some 
chain (5), hence where denotes the 

contribution of the chain (5). 

The principle of optimality states that the optimal 
chain (*) has the largest path integral at interval 
J — 1 as compared to all the other chains passing by the 
same chirplet in interval j. In particular, this means that 
ij-i = is larger that and 
Proof by contradiction: Let us assume that . 

We construct the chain (A) formed by (/3), the considered 
chirplet in interval j and the chain (<5). This CC is admis¬ 
sible. By construction, its path integral 
is larger than . Therefore, the chain (*) is not optimal 
which contradicts our hypothesis — QED. 

We apply this principle recursively starting from inter¬ 
val j = 0 and incrementing. For each chirplet interval 
and for all Nc chirplets of interval j, we keep only the 
CC maximizing the path integral up to this point and we 
discard the others. This procedure “prunes the combi¬ 
natorial tree” and avoids to consider useless candidates 
before going to the next interval. 

When the recursion reaches the last interval Nt — 1, 
we end up with a number W of CCs ending with a dif¬ 
ferent last chirplet and having the maximum path inte¬ 
gral among all chains with the same last chirplet. Fi¬ 
nally, within these “short-listed” candidates, we select 
the chain with the largest £ which is the global maximum. 


5. Numerical contingencies and computational cost 

To summarize, we started with the initial problem in 
Fq. inni) of finding the CC with the largest statistic. 
We rephrased this problem (using approximations) into 
a longest path problem in the TF plane. Here, path refers 
to the TF curve followed by the frequency of the CC, and 
the length is given by the integral of the WV of the data 
along the path. The maximization of the path length 


frequency 
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FIG. 4: Principle of optimality of DP 


over the set of CCs can be performed efficiently using 
DP. The resulting algorithm is tractable numerically as 
shown by the estimate obtained in the second part of this 
section. 

The definition of the CCs does not comprehend the 
fact that we only have access to discretized versions of 
the data and of their associated TF domain, denoted V 
in Sec. nm We begin this section by a discussion on 
these aspects. 

a. Discretization issues — On one hand, the defini¬ 
tion of the set of CCs relies on a TF grid sampling the 
continuous TF domain D. Theoretically, this grid can be 
refined arbitrarily. On the other hand, the search oper¬ 
ates effectively using the discretized version of I?, result¬ 
ing from the sampling associated to the WV. This fixes 
a maximum TF resolution which cannot be surpassed. 

It is useless to increase the resolution of the TF grid 
used for defining CCs beyond the one defined by WV. 
The WV divides the time axis into N intervals and the 
frequency axis into N bins^. Consequently, we have the 
following limitations, Nt < N and Nf < N. Further¬ 
more, in order to have time intervals (resp. frequency 
bins) of equal size , the TF grid parameters Nt (resp. 
Nf) must be divisors of N. 

All these requirements limit the choice of Nt and Nf. 
It may happen that the parameters of smallest tight CC 
grid are not suitable because of that. Note that in the 
case we consider, we are generally led to adopt the finest 
resolution for the frequency axis i.e., Nf = N. 

b. Estimate of the computational cost — We esti¬ 
mate of the computational cost by counting the float¬ 
ing point operations for all the primary subparts in the 
course of the procedure. The computation of the WV 
of the data involves N FFTs with time base 2Nf [^. 
such that the cost of this part is about 5 NNflog2Nf 
(assuming a standard implementation with RADIX-2). 

The number of operations required by DP is better 
estimated by grouping them by types, rather than by a 


^ It is possible to modify slightly the definition of WV in Eq. 1S3 
to get a finer sampling of the frequency axis and keep unitarity. 
We reserve this possibility for later investigations. 
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sequential assessment. The path integral ij in Eq. EH) 
is computed (with b additions) only once for each N^, 
chirplets of all Nt intervals, with a corresponding cost 
equal to NcN. 

For each of the Nc chirplets in each interval, the algo¬ 
rithm selects among the (at most) 2iV" -|-1 possibly con¬ 
nected paths. This procedure is repeated Nt — 1 times, 
and thus requires ~ NtNc{2N” + 1) operations. 

Knowing that the number of chirplets is Nc « (2A^r -I- 
1)A^/, the overall cost C thus scales with 

C cx bNNf\og^Nf+[N+{2N”+l)Nt]{2N'c+\)Nf, (63) 
which is a polynomial of the problem size. 

VI. APPLICATIONS 

In this section, the proposed method is evaluated with 
several numerical tests and compared with two other TF 
based algorithms for the detection of unmodelled chirps, 
namely the Signal Track Search (STS) m and Time- 
Frequency Clusters (TFC) 0 . The simulation code ® of 
these tests uses the implementation of these algorithms 
provided by |0. We Hrst give a brief presentation of 
STS and TFC. 


A. Existing algorithms 

1. Signal Track Search 

We have seen earlier that the TFR of a chirp signal 
can be essentially described in the TF plane as a regu¬ 
lar alignment of large values forming “ridges” along the 
instantaneous frequency evolution. The STS uses this ob¬ 
servation as a heuristic basis: detecting chirps amounts 
to hnding ridges in a TFR. 

In practise, the algorithm extracts the ridges from the 
WV distribution® of the data. Because of the presence 
of noise, image processing techniques are required to get 
a good ridge extraction. The authors chose an algorithm 
which is normally used for road extraction from aerial 
images. This algorithm is based on the fact that a ridge 
is a locus of points having a maximum curvature (as mea¬ 
sured by the second derivative) in the transverse direction 
and a small gradient along the longitudinal direction. A 
hysteresis thresholding procedure is applied over the sec¬ 
ond derivative of the WV (smoothed by a low pass filter) 


® Freely distributed scripts are available at 
http://www.obs-nice.fr/ecm for reproducing all the illus¬ 
trations presented here. 

9 In [13 , the authors use the standard definition of the disc rete 
WV originally proposed by Claasen-Mecklenbrauker (see for 
a definition and a detailed discussion). This definition differs 
from the one presented in Sec. |^,B_21 In particular, it does not 
satisfy unitarity. 


to detect TF points which suffice the above condition, 
and to grow iteratively chains of TF points from these 
ridge precursors. In [l2|, the ridge length (number of 
TF points in a ridge) is then employed as the detection 
statistic. However, we don’t use this definition here, but 
we rather consider the one given by the largest path in¬ 
tegral computed along the detected ridges. We observed 
that this variation outperforms the original definition of 
STS. 


2. TFClusters 

TFClusters is initially thought to detect short oscilla¬ 
tory transients (and not specifically chirps). The TFR of 
such transient is sparse i.e., the TF contents is essentially 
described by few components of large amplitude. The ba¬ 
sic idea of TFClusters is that, for reasonable SNR, the 
amplitude of the transient components is larger than the 
noise. 

This motivates the thresholding of the TFR of the 
data, given by the spectrogram (modulus square of short- 
time Fourier transform), to retain the TF points with 
the largest values. A clustering algorithm is used to 
group the selected points. “Signihcant” clusters are cho¬ 
sen whose cardinals are greater than a threshold. “In- 
signihcant” clusters are merged iteratively if they are 
sufficiently close to eventually form signihcant clusters. 
The statistic is then chosen to be the maximum sum of 
the TF powers over the clusters in the resulting list of 
signihcant clusters. 

3. Discussion 

It is important to stress a major difference between 
STS and TFClusters and the proposed method. By con¬ 
struction, the formers work well provided that the signal 
“stands above” the noise somewhere in the TF plane. If 
we dehne a local SNR in the TF plane (by computing at 
a given TF point, the squared difference of the mean val¬ 
ues of the TFR under the hypotheses Hq and Hi divided 
by its variance under Hq), then this is equivalent to say 
that the local SNR has to be large at least for some TF 
points. However, just like the standard matched hlter- 
ing, a detection with the best CC algorithm requires the 
global SNR (obtained by summing the local SNRs for all 
TF points) to be large. Clearly, this is a less stringent 
condition. 

TF path integration is a central ingredient of the best 
CC search. This idea is also used for other methods devel¬ 
oped for the detection of other GW sources, for instance, 
for inspiralling binaries, we can cite [^ and for the 
periodic GW sources^the Hough transform and the 
stack slide searches j^. 

Several distinctions must be stressed. First, the TF 
representation we use here (discrete WV) satisfies a spe¬ 
cific and crucial property, namely unitarity. This allows 
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us to link the final statistic to the quadratic matched fil¬ 
tering. TF representations based on short-time Fourier or 
wavelet bases used by the above methods are not unitary. 
Second, other methods require a precise model of the TF 
path (relying on the astrophysical source modelling) as 
opposed to our method. For the problem adressed here 
i.e., the detection of unmodeled chirps, we have shown 
that CCs can be treated as an effective finite template 
grid. We could then imagine to apply one of the above 
methods and integrate along the entire set of TF paths 
associated to CCs. This is however computationally im¬ 
possible because of the too large number of CCs, as al¬ 
ready discussed in Sec. IvU 

B. Newtonian chirps: illustrations and benchmark 

For the illustration of the best CC search, we use the 
Newtonian chirp signal introduced in Sec. IHTfI We 
recall that the frequency of such chirp is a power law 
given by Eq. Normally, the Newtonian chirp also 

includes a prescribed evolution of the chirp amplitude. 
However, for simplicity and better match with our initial 
model, we decide not to take this into account and keep 
the chirp envelope to a constant. 

The Newtonian chirp is completely defined by the total 
mass M of the binary (if we assume that the objects have 
equal masses) and its initial frequency /q. Fig. [^presents 
an example of a typical Newtonian chirp signal, where we 
set M = 7.3Mq and /o = 96 Hz. The chirp duration is 
T = 0.5 s. We fix the sampling frequency to fs = 2048 Hz 
(therefore, the number of samples is N = Tfs = 1024). 
A white Gaussian noise of unit variance is added to the 
signal. 

Within the GW literature, it is customary to define 
the SNR through matched filtering (assuming the initial 
phase is known a priori). We follow this definition which 
gives in the present case, 

N-l 

^^N/2. (64) 

k=0 

We note that, with this definition we have = 2£(s; (j)) 
(the factor of 2 accounts for the unknown initial phase). 

We choose to scale the chirp amplitude to a SNR p = 
20 . 

We apply the best CC search to this signal with the 
following search parameters. We arbitrarily fix the chirp¬ 
ing rate limits to be F = 8192 Hz/s and F = 917.5 
kHz/s^. These values are quite smaller than the ones 
expected at the LSCO (see Sec. IIV Ell but the time in¬ 
stant when theses limits are reached is close (few tenths 
of milliseconds before) to the LSCO. In Fig. 0 the 
time instants when the chirp (the solid line on the right 
panel) reaches the chirping rate limits (with dotted ver¬ 
tical lines) and when the binary system reaches LSCO 
(with dashed-dotted horizontal line) are indicated. We 
fix the frequency axis sampling to the finest accessible 


GW chirp in Gaussian noise, SNR=20 GW chirp [solid/green] best CC [dashed/red] 



time [s] time[s] 

FIG. 5: (Color online) Newtonian chirp in white Gaus¬ 
sian noise — left-. WV distribution of the signal. Only 
positive contributions are displayed (negative ones are set to 
zero) with a grey-scaled color map going from white (min¬ 
imum i.e, zero) to black (maximum), right: the best CC 
in dashed/red closely matches the actual instantaneous fre¬ 
quency in solid/green in the region where the regularity con¬ 
straints are satisfied. We indicate the instant when the chirp 
reaches the chirping rate limits with the dotted vertical lines 
and the frequency at LSCO with dashed-dotted horizontal 
line. 


resolution i.e., Nf = N = 1024. Similarly, we choose the 
smallest possible chirplet size with Nt = N/2. The rest 
of the parameters are derived from the regularity con¬ 
straints. In this respect, it is useful to calculate the adi- 
mensional characteristics of the problem i.e., N' = 2048 
and N” = 586.6. The resulting parameters are A/ = 9 
and A" = 3, which gives a maximum SNR loss p, « .28. 

We recall that the best CC search relies on the ap¬ 
proximation of the optimal statistic by a complex sum 
presented in Sec. Em The parameter p controls the 
relative precision of this approximation. From the re¬ 
sults of Sec. IV B H and the above general chirp specifica¬ 
tions, the approximation holds with a precision p = O.IA 
in a frequency bandwidth [/;, /s/2 — fi] with fi = 96 Hz 
which coincides (at least for the low frequencies, which 
are most important) with the frequency support of the 
present chirp. 

Fig. El presents the result of the best CC search with 
the above choice of parameters. The best CC closely 
matches the actual instantaneous frequency in the region 
where the regularity constraints are satisfied. 

An example is obviously not sufficient to evaluate the 
method thoroughly. Receiver operating characteristics 
(ROC) gives a systematic assessment of the performance. 
The ROC of a given statistic I is the diagram giving the 
detection probability Pd{lo) = ]P(^ > ^o|Ai) versus the 
false alarm probability Pfa{lo) = P(^ > ^o|Ao) at a given 
SNR and for all thresholds Iq. 

For this exercise, due to computing limitations, we pre¬ 
fer short signals with a small number of samples A. We 
choose a Newtonian chirp with total mass M = IIMq 
and initial frequency /o = 96 Hz which has a short du¬ 
ration T ~ 250 ms. Choosing the sampling frequency 
fs = 1024 Hz, we have A = 256 samples. White Cans- 
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sian noise is added to the signal and the amplitude is 
scaled such that the SNR is p = 10 . 

We fix the chirping rate limits to F = 8.192 kHz/s 
and F = 1.05 MHz/s^. Like the above example, these 
limits are reached at a time instant close to the LSCO. 
We choose the finest TF grid parameters Nt = 128 and 
Nf = 256 , and the regularity parameters N^. = 9 , iV" = 
4. The resulting CC grid is tight with p. « 0 . 4 . 

Concerning STS^° and TFClusters^^, we set their free 
parameters empirically using the recommendations avail¬ 
able in the references, without a precise fine-tuning. Fig. 
El displays a single trial and Fig. [7| presents the ROCs of 
the three methods presented previously. We see that the 
best CC search outperforms the two others as expected. 

Here, we wish to add few remarks regarding the com¬ 
parison between the best CC search and STS. The im¬ 
provement in the ROC of the best CC with respect to 
STS has two origins. First, the use of a unitary discrete 
WV instead of the standard WV helps in increasing the 
detection probability by few percent. The unitarity pre¬ 
serves the power in TF plane and hence improves the 
efficiency. Second, the major part of the improvement 
comes from the TF pattern search procedure. As ex¬ 
plained in Sec. IVT^ the use of a global search criterion 
instead of a local one is a crucial ingredient. 

It is interesting to compare these ROCs with what 
could optimally achieve an imaginary observer which 
knows in advance the targeted chirp. Since this clair¬ 
voyant observer knows the chirp phase exactly, he can 
apply the optimal statistic i.e., the quadrature matched 
filter obtained in Appendix IXI The ROCs of the quadra¬ 
ture matched filter can be obtained analytically (under 
Gaussian noise hypotheses). The false alarm and detec¬ 
tion probabilities are given respectively by [^ : 

Pfaik) = exp(-Zo), ( 65 ) 

P,{lo) = 1 - eM-plm E in + 1) (66) 

n\ 

n=0 

where Iy{x) = l/r(a;) du is the incomplete 

Gamma function. 

This ROC depends only on one parameter, namely the 
SNR Pc- The ROC curve of clairvoyant statistic with 
Pc = P provides an absolute upper bound on the detection 
probability. Obviously, having in hand all the informa¬ 
tion makes a very large difference with respect to the case 


Following the notations of , the size of the Gaussian kernel 
of the pre-smoothing filter is fixed to cr = 2. The low and high 
thresholds of the hysteresis are set to 3.3/pixel^ and 10/pixel^ 
resp. 

The TFR is given by the short-time Fourier transform computed 
over non-overlapping blocks of 16 samples (i.e., intervals of ^ 7.8 
ms). The frequency axis is tiled into 32 bin s (i .e., a resolution of 
32 Hz). We use the nominal values given in Il4ll for the rest of the 
parameters namely p = 0.1, fj = 5, 5 = [0, 0, 0, 0, 0, 0, 2, 3, 4,4] 
and a = 0.25. 


GW chirp in Gaussian noise, SNR=10 
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FIG. 6: (Color online) Newtonian chirp in white Gaus¬ 
sian noise — left: WV distribution of the signal (displayed 
similarly as in Fig. El- right: actual chirp frequency in 
soiid/green and best CC in dashed/red. We indicate where 
the chirp reaches the chirping rate limits with dotted vertical 
lines and the frequency at LSCO with dashed-dotted horizon¬ 
tal line. 


where we only know that the incoming GW is a smooth 
chirp. The detection probability of the clairvoyant statis¬ 
tic is very close to I over the entire range of value chosen 
for the false alarm rate. This is why we don’t show this 
curve. It is more interesting to compare the performances 
of the various methods with the ones of the clairvoyant 
observer for SNRs pc < p. More precisely, we adjust pc 
in such a way that the resulting curve matches reason¬ 
ably well the ROC of the best CC search in the region 
of interest i.e., for false alarm probabilities in the range 
10“^ to 10“^. Since the SNR is inversely proportional 
to the distance of the GW source, the ratio of the actual 
SNR to the best-fit value pjpc gives the reduction factor 
of the sight distance with respect to the ideal (and non 
accessible) situation where we have at our disposal all 
the information about the chirp we want to detect. We 
include the fitted clairvoyant ROC in Fig. 0 The ratio 
in the sight distance can be estimated ~ 10/6.15 « 1.6. 


C. Random CCs and robustness 

While benchmarks based on Newtonian chirps are sat¬ 
isfactory for a comparison of several detection methods 
in a nominal situation, they don’t provide a test for the 
robustness i.e., a measurement of their ability to detect 
reliably a large class of different chirps. 

In this section, we present ROC curves computed us¬ 
ing random CCs. Random CCs are generated by chaining 
chirplets randomly chosen in a range specified by regular¬ 
ity constraints. Therefore, the frequency of a random CC 
follows a kind of random walk in the TF plane. We gen¬ 
erate a new random CC for each trial made to estimate 
the detection probability. 

The detection of random CCs is obviously much more 
difficult than the detection of a single chirp. It is an 
effective test of the method robustness. No classical ap¬ 
proaches (e.g., based on banks of quadrature matched 
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receiver operator char. (N=256, SNR=10.0) 



FIG. 7: (Color online) Newtonian chirp in white Gaus¬ 
sian noise — Comparison of ROCs of the best CC search 
(dashed/blue, with error bars in solid/red) with STS (dot¬ 
ted/magenta) and TFC (dashed-dotted/cyan). The compn- 
tation of each ROCs is perfomed over 2 x 10® trials (half for 
the false alarm probability and half for the detection probabil¬ 
ity). The diagram also includes the ROC of the clairvoyant 
quadrature matched filter (bold dashed/green) shown here 
with the SNR pc = 6.15 adjusted to reasonably fit the ROC 
of the best CC search. 


filters as for inspiralling binary chirps) can be applied 
successfully in this case. 

We assume the same general characteristics of the 
Newtonian chirp used in the first example in the previous 
section, namely T = 0.5 s, /s = 2048 Hz, thus N = 1024 
samples, F = 8192 Hz/s and F = 917.5 kHz/s^. We 
already computed satisfactory search parameters for this 
set-up. Therefore, they remain unchanged (W = 512, 
Nf = 1024, TV/ = 9 and iV" = 3). The random CCs 
are generated on the same basis, but with a time interval 
slightly larger, the regularity parameters being increased 
accordingly i.e., W = 64, TV/ = 1024, TV/ = 65 and 
TV/' = 57. We use an additive white Gaussian noise. 

Fig. E presents an example of such signal (with SNR 
p = 20) and the result of the application of the best 
CC search. Fig. El displays the ROC curve of the best 
CC search (with SNR p = 12) along with the one of 
the clairvoyant quadrature matched filter adjusted to an 
adequate SNR. We estimate a loss in the sight distance 
with respect to the clairvoyant case to be a factor of 
~ 2.6. Best CC search “sees” to distances comparable to 
(in the sense, with a reduction factor less than one order 
of magnitude) what classical methods achieve in other 
GW detection problems. 

The computational cost of this search as estimated by 
Eq. (IHHll is about 142 millions of floating points opera¬ 
tions for one block of duration T — 0.5s. Assuming 10% 
overlap between successive blocks, real-time processing 
can be achieved with a computing power of 2.8 Gflops 
which is less than what a single standard workstation 
can handle today. 
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FIG. 8: (Color online) Random CC in white Gaussian 
noise — In these plots, we arbitrarily set fs = 1. top - left: 
example of a random CC in white Gaussian noise, top - right: 
noise free random CC. — bottom - left: WV distribution 
of the signal (displayed similarly as in Fig. EJ. bottom - 
right: actual chirp frequency in solid/green and best CC in 
dashed/red. It is worthwhile to note that, although the best 
CC can lose track for some time because of noise fluctuations, 
it is able to recover the exact TF path. 
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FIG. 9: (Color online) Random CC in white Gaus¬ 
sian noise — This diagram displays the ROC of the best 
CC search (dashed/blue, obtained from 2 x 10® trials) com¬ 
pared with the analytical ROC of the clairvoyant matched 
filter (bold dashed/green) with the SNR pc = 4.55, adjusted 
to produce a reasonable fit. 


VII. CONCLUDING REMARKS 

Smooth chirps define a general model of “nearly phys¬ 
ical” GW chirps. Chirplet chains - chains of linear 
chirplets - allow the design of tight template grids for 
the detection of smooth chirps. The optimal detection 
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requires these grids to be searched thoroughly to find the 
template which best matches with the data. Although 
the shear large number of templates prevents the use 
of an exhaustive search, near-optimal detection can be 
performed with the time-frequency based procedure pre¬ 
sented here. Its originality lies in the clear link estab¬ 
lished between the optimal statistic and the proposed 
search algorithm as opposed to with other approaches. 
In particular, it justifies the choice of a specific time- 
frequency representation (the unitary discrete WV) and 
pattern search algorithm (dynamic programming). We 
have evaluated that best CC search is computationally 
tractable for detection of typical GW chirps. 

It is important to emphasize several features which 
makes the proposed method attractive in practise. First, 
the free parameters (the chirp duration T and the chirp¬ 
ing rate limits F and F) are few and directly related to 
physical characteristics. Second, the principle “77e who 
can do more can do less" applies here: smooth chirps is 
a very general class of chirps. This model, and thus the 
search algorithm can be easily modified and adapted to 
incorporate additional astrophysical information. For in¬ 
stance, it is easy to search only chirps with an increasing 
(or decreasing) frequency. One may also want a more 
stringent constraint on the chirping rate at low frequen¬ 
cies than at high frequencies. The inclusion in the algo¬ 
rithm of a dependency of the chirping rate limit upon the 
frequency is straightforward. This leaves the possibility 
of a compromise between efficiency (since the restriction 
of the set of admissible waveforms due to additional con¬ 
straints reduces the false alarm rate) and robustness, de¬ 
pending on the quantity and reliability of the information 
available on a specific GW source. Third and finally, it 
is simple to restrict the search to chirps starting and/or 
finishing at given time-frequency location. This feature 
could be used for partially known chirps whose waveforms 
is known only on a part of the total duration. Those sig¬ 
nals could be detected with a hybrid approach combining 
a standard matched filtering where the waveform model 
is available, and best chirplet chain search for the rest. 


Acknowledgments 


recall that 

1 

A{x; {A, ipo, to, />(•)}) = ^ X! 

\fc=o / 

where A/" = norm of Sk = cos{(j)k+‘Po)- To 

keep the notations simple, we don’t mention all parame¬ 
ters explicitly and set A{x; (po) = A(x; {A, ipQ, to, </(•)})■ 

In the literature concerning the detection of inspi- 
ralling binaries of compact objects min , this maxi¬ 
mization is usually performed assuming that Af is inde¬ 
pendent of ipo- This assumption is correct when the two 
quadratures cos (f>k and sinc/fe, viewed as vectors of 
are orthonormal (i.e., orthogonal and of same norms). In 
this case, we have Uc = Ug = N/2 and = 0 where 
ric, ns and Ux are the norms and cross-products of the 
quadratures as defined in Eqs. O- Inserting this into 

N = nc cos^ (/?o ~ nx sin(2(po) + ng sin^ (/?o, (A-2) 

we conclude that Af = N/2 is a constant. 

However, for general phase evolution, the quadrature 
waveforms are not necessarily orthonormal. This is ap¬ 
proximately true when the chirp oscillates sufficiently 
rapidly during many cycles (e.g., for inspiralling binaries 
of a small mass). Since we are considering chirps with an 
arbitrary phase and of relatively short duration, such as¬ 
sumption is not realistic and we opt for the general case 
keeping the dependency of Af upon po. 

Expanding Sk in terms of two quadratures and rewrit¬ 
ing Eq. m, we get 

(a:cCOS(/?o - a;sSin(po)^ 

■^(^7 V^o) , 2 • /O ^ I * 2 ’ 

2[nc cos^ (fo — Tlx sm{2ipo)ris sm (po) 

(A3) 

where Xc and Xg are the cross-correlation of the data with 
COS (j)k and sin (fk as defined in Eq. ®. 

To proceed with the maximization, we first examine 
the special case where the quadratic waveforms are lin¬ 
early dependent i.e., cosc/k oc sinc/fe for all k. This im¬ 
plies that we are in the degenerate case where ipk = Po is 
constant. Introducing the two angles p = arg(a;c -I- iXg) 
and r] = aTg{y/nf + we can rewrite Eq. 1A3II as 
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A(x;(po) 


Ix1 + x1 cos'^{p + po) 
2nc + Ug cos2(?7 -I- po) 


(A4) 


The proportionality of the quadrature waveforms im¬ 
plies that ^ynfxg ± ^/nfxc = 0 which gives s\n.{ri — p) = 0, 
and hence rj = p + ttZ. We conclude that A{x;po) re¬ 
mains constant for all po and is equal to the statistic 
given by 


APPENDIX A: MAXIMIZING LLR A( ) OVER 
THE INITIAL PHASE po 

In this appendix, we maximize the statistic 
A(a;; {A, (po, to,'/(•)}) over the initial phase po. We 


£(x;to,</) = ^^^. (A5) 

In the non degenerate case, we compute the derivative 
of the statistic as given in Eq. (IA3I1 w.r.t. pq. Its numer¬ 
ator turns out to be a second order polynomial of tan pQ. 
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The root associated to the local maximum is 

— 1 f 


(pQ = tan 


TlxXa TlaXc 


(A 6 ) 


at ^ ^ \ nsxl-2n:,XcXs+ncxl 

£{x;to,(l)) = Aix;ipo) = - ^ -, (A7) 


which gives the ML estimator of the initial phase. 
Inserting this expression in Eq. (IA3II yields 


where O = ncUs — > 0 . 

We can re-express this statistic as 

\ 2 


i{x-,to,(l)) = i 


/N-l \ 

^ ^ Xk^k 
\k=0 / 


N-1 \ 

^ ^ Xk^k 
k=0 / 


(A 8 ) 


where Ck and Sk are the orthonormalized counterparts of 
the waveforms in quadrature cos (ftk and sin (j)k obtained 
from the Gram-Schmidt procedure as given below 


cos 4>k 



Sk 


ric sin (j)k - cos 4>k 
VncO 


(A9) 


and referred to as templates of (j). 

In practise, this orthonormalization is indeed per¬ 
formed for the detection of inspiralling binaries (see , 
p. 3046) and is justified with heuristic arguments. The 
derivation shows that it results directly from the maxi¬ 
mization of the LLR. 


APPENDIX B: TAYLOR APPROXIMATION OF 
THE DISTANCE BETWEEN CHIRPS 

In this appendix, we detail the approximation of the 
statistic £{s] </>*) with Sk = Acos((j)k + <po) and assuming 
that the template phase (j>* is close to the phase (p of the 
signal s present in the data. We start from the following 
Taylor expansion of i{s; p*) for small Ak = 4>l — (pk 

N-l 

e{s-,cp*)=Ks-A) + Y. 9ke\<p^=^Ak 

k^O 

N-l 

+ 2 ^ dkA<p*=rpNAi + ..., (Bl) 

k,l—0 

where the partial derivatives dk = djQfpX, and = 
ldp\dp\ are taken with respect to the samples of the 
template phase (p*. Next, we examine this expansion 
term by term and obtain analytical expressions as a func¬ 
tion of the phase samples {pk} and {PX\- 


1. First derivative: local extremum 


We get the following general expressions of the deriva¬ 
tive of the numerator 

Pkkl — dk^sXf, Tig2Xc.dk^c ‘2{^3kklxXcXs T klxdkXcXg 

TlxXcdkXs) ^k'klcXg ‘kl(p2XsdkXs^ 

and of the denominator 

dkd = 2{dkncns + ncOkUs - 2nxdknx)- (B3) 

We insert Sk = Acos{pk + V^o) and work out each of 
their component term. At the match (when p* = p), we 
get 


(/) — Sin2pk , 3kklx\(l)*—(j) — cos2pk , (ld4) 

Pkklclcp*—(j) — Syil2pk — 3kk^s\(p*—(p ^ 

dkXs\ 4 ,^=p, = AcospkCosipk + ipo) ■ (B6) 

dkXc\p,*=^ = -Asini^fc cos{pk + <po), (B7) 

Combining all the above expressions, the derivative 
dkiT, can be factorized, yielding 

3kk^\(p*—<p — "^(^j P^3kd\^*—(p , (B8) 

where dkd\^*=^ = 2{nc — ns)sm2pk — iux cos2pk and 
£{s]p) = A‘^{nc cos^ (fo — Ux sm2(po + Us sin^ (fo)/2. In 
conclusion, the first derivative 9fc^|0*=(^ = 0 vanishes at 
p* = p which is thus a local extremum. 

Using the parameters e = 2nx/N and 6 = (ric — ns)/N 
as defined later in Sec. IV Bl (and also discussed in Ap¬ 
pendix P, the statistic and the denominator at the 
match can be expressed as functions of 5 and e as 

A^N 

£{s]p) = — - —(1-I-<5cos2(/Jo — esin2y)o) (B9) 
dU*=0 = —(1-^^-e^)- (BIO) 


2. Second derivative and distance 


We show in the previous subsection that the first 
derivative at the match dk£\ 4 ,‘= 4 , vanishes. Consequently, 
the second derivative at the match can be expressed sim¬ 
ply in terms of the second derivatives of the numerator 
and denominator at the match, namely 


9ki^\<l>*=ij> — 


dhn- £{s-,p*)dlid 


(Bll) 




We obtain the following general expressions for the sec¬ 
ond derivatives of the denominator 


From Eq. ®, we write the statistic £ as the ratio £ = 
n/d. The numerator is n = ngx'j. — 2nxXcXs + ricXg and 
the denominator is d = 2{ncns — n'^)- We thus have 
dk£ = {dkn - £dkd)/d. 


^kld — ‘^\Pkl^ck^s T Pkk^cPlk^s T Pl'kl'cPk'kl's T klf^dj^iTlg 

- 2{dknxdinx + nxdhnx)]. (B12) 

and of the numerator 













19 


Here, we adopt the precedence rule df^ab = {df,a)b. 


^kl^s^c ‘^i,^k^s^c^l^c ~t“ ^l^s^c^k^c ~t“ klgdj^X qDiXq -\- TlgXQOf^iXo^ 

s ~\~ dj^TlxOlXcXs df^TlxXcdlXg “t“ dlTlxdk^cXs “t“ dlTlxXc^kXs 

-\- Tlxdl^XcdlXs klxdlXcdkXs TlxXgdj^iXc “t“ TlxXcd^iXg') 

^kl^c^ s ‘^(^^kklc^s^l^s “t“ dlTlcXgdf^Xs kl^dj^XsOlXs “t“ TIqX gOf^^X s') i 


Similarly to the first derivative, we insert the expres¬ 
sion of the signal Sk = Acos{(j)k -I- ^o) and evaluate each 
of the component terms of the above expressions. We 
have to distinguish two cases i.e., the non-diagonal cross 
terms of the Hessian matrix when k ^ I and the diagonal 
ones when k = 1. 


derivatives are zeros (namely d^iric = dhux = dhug = 0 
and dhxc = d^Xg = 0). We get 


=-4cos2(0fe - (/>;), (B14) 


a. Cross terms, k A I 


When k ^ I, the above Eqs. (ITO and iirm are 
significantly simplified because all the second order cross and combined with Eq. (15^ . 

I 


['9fezn - i{s; (l)*)dlid],j,*=^ = [1 -h cos 2(<))k - (/>i) - cos 2(<))k + ipo) - cos 2(^; -h <f>o) 

+ e (sin2(^fc + + ‘Po) — sin2(/)fe — sin2^j — sin2(/5o) 

-l-i5 (cos 2((/)fe -I- -I- (po) — cos 2(j)k — cos 2(j)i + cos 2(/5o)] (B15) 


In Appendix m we discuss the range of values taken 
by e and 5 depending on the phase (f>. We show that 
these parameters are small e,S <St. 1 if the phase </) is a 
CC whose frequency does not come close to DC nor the 
Nyquist frequency. We assume that this remains true in 
the more general case, when cj) is the phase of a smooth 
chirp. We retain the leading term (of order 0 in e and S) 
and get the following approximation 


dlA^>=4> = Xki = ^ ((1 - Cfe)(l - Cl) + Skh) (B16) 


where Ck = cos2{(j)k + </^o) and Sk = sin2{(j)k + <Po)- 


b. Auto terms, k = I 

We consider the case where k = 1. Now, the second 
order derivatives do not vanish. In fact, we have 

dkric\(p*—(p — 2 cos 2<()/j; dkrix\cp*—(p — 2s\Yi2(f)k 

(BI7) 

dlng\^»=^ = 2cos2(j)k (BIS) 

dlxc\^*= 4 , = -A/2{cos{2(j>k + ipo) + cos (po) (BI9) 

9fca;s|0*=0 = -A/2{sm{2(j3k -k (po) -sint/Jo). (B20) 

The consequence is an additional term Du to the sec¬ 
ond order derivative of the statistic = Xu + 

Du- With a direct calculation, we obtain its exact ex¬ 
pression (no approximation needed): 

y\2 

Du =+ Ck)6u- (B21) 

where the Kronecker symbol is Su = 0 for k ^ I and 1 
for k = 1. 
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3. Approximated distance 

From Eqs. fFTTlTl and assuming that |e|, |<5| 

1, we have £{s;(j)) ~ A^N/A. The distance defined in 
Eq. (11311 can thus be written as 

JV-l / N-l \ 2 

^ E (1 - E (1 + 

k=0 V ^ k=0 ) 

1^-1 

^ E! ■Sfc^fe I ■ (B22) 

^ k=0 ) 

Considering that and are slowly varying with 
respect to c/c and Sfe, we argue that, similarly to what is 
discussed in AppendixIHJ the positive and negative terms 
compensate when making the following sums ’^^CkA.k, 
SfcAfe and X^fc^fcA^. We neglect the small residual, 
which leads to the final approximation of the distance in 
Eq. itni) . 


the identification of the maximum) and concentrate on 
“cross-terms” (i.e., terms in UjUk)- A direct calculation 
leads to 

N-2 N-l 

V{r) = Va + E E 

i=l k=j+l 

where cjk = 2j{N — k)/N^ and Va is the contribution 
due to the auto-terms. 

Since all Cjk > 0, the maximum of V is reached when 
all Uk have the same signs, that is when Uk are all 
identically +U or —U. Therefore, the empirical vari¬ 
ance is maximum when = ±kU and in this case 
V{r) = U^{N^ - 1)/12. 

We recall that the distance between the smooth chirps 
is well approximated by the empirical variance of the 
phase discrepancy [see Eaini. We apply this result to 
the original maximization problem by setting rfc=Afe and 
U=2TrAfts as given in Eq. (ITTIl . 

APPENDIX D: BOUNDING S AND e OF A CC 


APPENDIX C: CONSTRAINED MAXIMIZATION 
OF THE DISTANCE 

We rewrite the constrained maximization problem de¬ 
scribed in Sec. II V C 31 of the distance in Eq. m under 
the constraint in Eq. with simpler notations. We 
relate them to the initial problem at the end of this Ap¬ 
pendix. 

Let {rk} a series of N real numbers. We want to max¬ 
imize the empirical variance V (r) expressed by 

iV-l / iv-i \ 2 

no = ^E^E ^E^M ’ (Cl) 

fc=0 \ k=0 / 

under the constraint that the increments Uk = rk — rk-i 
are absolutely bounded by some constant U > 0 i.e., 
|ufc| <U for fc > 0. 

The empirical variance V(r) is invariant by the addi¬ 
tion of an arbitrary constant C: let = yk + C, for all k, 
then V{r) = V{y). We can thus assume with no loss of 
generality that tq = 0 (i.e., choose C = — j/o)- Therefore, 
we have rt = for fc > 0. 

We want to maximize the convex function V in the set 
of feasible solution described by {r^} which is a polyhe¬ 
dron of Erom a classical theorem of convex analysis 
(see lU, P- 187), we conclude that V reaches its max¬ 
imum at one of the extreme points of this polyhedron. 
The extreme points are the points where the increments 
are either Uk = +U or Uk = —U. There are 2^“^ extreme 
points and we need to identify the one which maximizes 
the convex function. 

Let us rewrite the empirical variance V (r) as a function 
of Uk- We leave the “auto-terms” aside (for all extreme 
points, the auto-terms are equal to independently of 
the sign of Uk- Their contribution is thus unimportant for 


The simplification of the statistic in Sec. lVBl is closely 
related to the orthogonality and length difference of 
the vectors c = {cos^fc,fc = 0, ...,A — 1} and s = 
{sin(j)k, k = 0,..., N — 1} of . 

Noting that their norms and scalar-product are respec¬ 
tively given by n-c = (c, c), ris = (s,s) and = (c,s) as 
defined in Eq. o, the departure from “orthonormality” 
of c and s can be quantified by the two parameters 


ric + ns nc + Us 

The parameter S measures the relative difference of the 
lengths of c and s while e is related to the cosine of the 
angle between the two vectors. 

When the vectors c and s are orthonormal i.e., orthog¬ 
onal and of same lengths, both S and e are zero. By 
continuity, for nearly orthonormal vectors, 6 and e are 
then expected to be small. Intuitively, this should be 
true for vectors with oscillating components like c and s. 
Indeed, e and S can be rewritten in the form of oscillating 
sums, namely 

^ Af-l ^ N-l 

S = —cos2(l)k e = — '^sin2<j)k. (D 2 ) 

fe=o fe=o 

The positive and negative contributions cancel in the 
summation, and thus leaves a small residual. In this ap¬ 
pendix, we go beyond this intuitive rationale when the 
phase ^ is a CC as defined in Eq. m and give a system¬ 
atic investigation of the maximum value taken by 6 and 
e. 

Eq. dni motivates us to combine 6 and e is the fol¬ 
lowing complex sum S 

I w-i 

S = S + ie = — expi 2 (^fe. (D3) 

k=0 
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Bounding the modulus of S is equivalent to bounding S 
and e. Analytic number theory provides a large number 
of results concerning exponential sums like S, for improv¬ 
ing upon the trivial bound jS”! < 1. We use one of these, 
namely the Kuzmin-Landau lemma, see |32| p. 7. We 
present a proof of this lemma pertaining to the present 
case where the phase (j> is a. CC. 

The proof can be summarized as follows. A change of 
variables is introduced which allows us to put a bound 
on the modulus of S' by a sum of the finite difference 
of complex variables. These new variables appear to be 
collinear in the complex plane. The sum of the modulus 
of their difference is thus equal to the distance between 
the extremes. The final bound on |S| is then obtained by 
combining this property with the explicit expression of 
the phase of the CC, provided a constraint on the lower 
and higher frequencies reached by the CC. 

Let us define for 1 < fc < — 1, the following variables 

(D4) 

We perform the above change of variables in the sum 
S using the relation 

exp{i24>k) = [exp(z2^fc) - exp{i2(j)k+i)]Ck+i , (D5) 

and we get 


N-2 

NS = Cl exp(*2(/)o) + ^ (Cfe-i-i - Cfe) exp{i2(l>k)+ 

fe=i 


(1 - CAr-i)exp(z2C)Ar-i). (D6) 


By taking the modulus on both side, we obtain the 
following bound. 


Nt-l U+l)b-l 

A^|S| < |Ci| + |l-Cw-i| + E E ICfe-i-i ~ Cfcl+ 

j—0 k—jb-\-l 
Nt-2 

E! ICj'ft-i-i ~ CiftI- (D7) 

where we split the sum in Eq. inn) into smaller ones 
calculated over chirplet intervals i.e., tj < tsk < tj+i 
or equivalently jb < k < {j + l)b — 1 with b = St/ts, 
the number of samples in a chirplet interval. In the last 
sum, we separate the terms corresponding to the tran¬ 
sition between two consecutive chirplets from the terms 
corresponding to the individual chirplets. 

We now obtain a bound on each term of RHS of 
Eq. (inTll . starting with the first sum. Eq. fni) can be 
rewritten as 


Cfe = ^[1 + *cot((ifc/2)]. (D 8 ) 

The variables Cfc are all located on the line Ift(C) = 1 / 2 . 


Within a chirplet interval, i.e. if tj < ktg < tj+i, the 
phase difference is a linear function of k given by 

dk = 2Trts [(2 - r)f^. + , (D9) 

where r = {2tj^k — ts)/St- 

We assume that the node frequencies of the CC are 
constrained in the following bandwidth : 

fl < fm, < fs/2 - fl, (DIO) 

where fi = fsC/2 and 0 < c < 1/2. In other words, the 
CC cannot approach arbitrarily close to neither DC nor 
the Nyquist frequency. 

Since 0 < r < 2, we have ‘iirtsfmj S dk < if 

fmj < frrij+i (and the opposite in the other case) which 
implies that 


0 < 27rc < dk < 27r(l — c) < 27r, (Dll) 

for all fc, hence —oo < ^(Cfc) < +oo. 

If fnij < frrij+i (resp. fmj > /m^+J, the phase dif¬ 
ference dk and hence 3(Cfc), increases (resp. decreases) 
monotonically with k. 

Since their imaginary parts are finite and monotonic, 
the variables Qk are associated to consecutive points on 
the line 3fi(C) = 1/2 of the complex plane. The sum of 
the lengths of the segments linking two nearby points is 
equal to the length between the extremes, thus 

(j+i)b-i 

E ICfe-i-i - Cfcl = ICa-i-i)& - Oft+il ■ (D12) 

k=jb+l 

Applying the mean value theorem to the func¬ 
tion g{x) = cot(a;/2)/2, whose derivative is g(x) = 
l/(4sin^(a;/2)) and using the constraint in Eq. iiTrrni . 
we obtain the following bound 

lA /■ \ ^ \d{j+l)b - djh+i\ 

ICO'+l)b G'b+ll ^ A ’ 2/ \ ■ 

4 Sin (ttc) 

We carry on by bounding the numerator 

|d(j+i)b - djb+i I = ^'r^ts\fm.j+i - fmj I < (D14) 

and denominator with 2c < sin(7rc) (this is valid for 0 < 
c < 1/2) and by summing over all j to finally obtain the 
bound on first summation term in Eq. 


Nt-l 

E 


U+i)b-l 

E iCfc-i-i - Ck\ < 

k=jb 


irN'^Nt 
ANc^ ■ 


(D15) 


The second summation coming from the boundary 
points of the chirplet intervals can be bounded in a sim¬ 
ilar way, considering that 


ICjb+l Ob I ^ 


|dj5-t-l bljb\ 
4 sin^( 7 rc) 


(D16) 
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and combining with 

\djb+i — djb\ = 27rtg|/mj_,.i —/mj-i |/<5t ^ 


2 ,. 2 i(D17) 


we get the result 

Nt-2 


iCjfc+i Oft I ^ 


i=i 


TTtsKNt 

AN c^Sf 


(D18) 


Finally, from Eq. irna . we have the following inequal¬ 
ities 


This bound is obtained from a worst case estimate. 
Generally, 5 and e are smaller than this value. With the 
choice of a small c, a more realistic estimate rather than 
a strict bound can be obtained replacing the inequality 
2c < sin( 7 rc) by the first order Taylor approximation ttc ^ 
sin(7rc) in the proof above, yielding the following estimate 


KNt 

7rg2jY2 


(D22) 


2| sin(dfe/2)| 2 sin( 7 rc) 4c’ ^ 

which, when applied with k = \ and k = N — 1, set 
an upper limit to the remaining terms in the RHS of 
Eq. 1D7II . noting that |1 — Cn-i\ = |CAf- i|- 

Combining this result with Eqs mi511 and (ID18II . we 
get 


Summarizing, we obtained an upper bound 77 on IS”! by 
restricting the frequency of the CC in a bandwidth de¬ 
fined by c. We rather use the reciprocal i.e., we get the 
limits of the frequency bandwidth from an acceptable 
value for 77 . If we assume that ~ A{N'/Nt){Nf /{2N)) 
as given by the regularity condition, the frequency band¬ 
width is [/i, /s /2 - /;] with 


1^1 < 


1 

'Wc ^ 


irN'^Nt 

4c2Af2 


(1 + 1 /&). 


(D20) 


The number of samples in a chirplet interval being an 
integer 6 > 1, and selecting the dominating contribution, 
we conclude that IS"! < 77 with 


_ TTiV/iVt 
^ “ 2 c27V2 ■ 


(D21) 


fi 


fsC 

2 


2 .hVN'6f 




(D23) 


where the leading constant is obtained from ^J2^Jti 
2.5. We use this result in Sec. mu 
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